#include <iostream>
#include <opencv2/core.hpp>
#include <opencv2/imgcodecs.hpp>
#include <opencv2/imgproc.hpp>
#include <opencv2/highgui.hpp>

#define EPSILON 1E-5

static double _maximum(double num1, double num2, double num3) {
    return std::max(std::max(num1, num2), num3);
}

static bool _almost_equal(double num1, double num2) {
    return (std::abs(num1-num2) <= (EPSILON * _maximum(1.0, std::abs(num1), std::abs(num2))));
}

static bool _greater_or_equal(double num1, double num2) {
    return ((num1 > num2) || _almost_equal(num1, num2));
}

static unsigned int _successor(unsigned int idx, unsigned int numOfPoints) {
    if (idx+1 == numOfPoints) {
        return 0;
    }
    return idx+1;
}

static unsigned int _predecessor(unsigned int idx, unsigned int numOfPoints) {
    if (idx == 0) {
        return numOfPoints-1;
    }
    return idx-1;
}

static double _distance_from_point_to_line(const cv::Point2f& pt, const cv::Point2f& linePtA, const cv::Point2f& linePtB) {
    double term1 = linePtB.x - linePtA.x;
    double term2 = linePtA.y - pt.y;
    double term3 = linePtA.x - pt.x;
    double term4 = linePtB.y - linePtA.y;
    double nominator = std::abs(term1*term2 - term3*term4);
    double denominator = std::sqrt(term1*term1 + term4*term4);
    if (denominator != 0) {
        return nominator / denominator;
    }
    return 0;
}

static double _height(const cv::Point2f& pt, const std::vector<cv::Point2f>& polygon, unsigned int c) {
    unsigned int numOfPoints = polygon.size();
    cv::Point2f ptC = polygon[c];
    cv::Point2f ptCPre = polygon[_predecessor(c, numOfPoints)];
    return _distance_from_point_to_line(pt, ptC, ptCPre);
}

static double _height(unsigned int ptIdx, const std::vector<cv::Point2f>& polygon, unsigned int c) {
    //int numOfPoints = polygon.size();
    cv::Point2f pt = polygon[ptIdx];
    return _height(pt, polygon, c);
}

static void _advance(unsigned int& idx, unsigned int numOfPoints) {
    idx = _successor(idx, numOfPoints);
}

static void _advance_b_to_right_chain(const std::vector<cv::Point2f>& polygon, unsigned int& b, unsigned int c) {
    unsigned int numOfPoints = polygon.size();
    unsigned int succ = _successor(b, numOfPoints);
    double succBHeight = _height(succ, polygon, c);
    double bHeight = _height(b, polygon, c);
    bool greaterEqual = _greater_or_equal(succBHeight, bHeight);
    while (greaterEqual) {
        _advance(b, numOfPoints);
        succ = _successor(b, numOfPoints);
        succBHeight = _height(succ, polygon, c);
        bHeight = _height(b, polygon, c);
        greaterEqual = _greater_or_equal(succBHeight, bHeight);
    }
}

//////////////////////////////////////////////////////////////////////////////////////////////////////////////

static bool _are_equal_points(const cv::Point2f& pt1, const cv::Point2f& pt2) {
    return (_almost_equal(pt1.x, pt2.x) && _almost_equal(pt1.y, pt2.y));
}

static void _line_equation_by_points(const cv::Point2f& pt1, const cv::Point2f& pt2, double& a, double& b, double& c) {
    CV_Assert(!_are_equal_points(pt1, pt2));
    a = pt2.y - pt1.y;
    b = pt1.x - pt2.x;
    c = -a*pt1.x - b*pt1.y;
}

static std::vector<double> _line_equation_parameters(const cv::Point2f& pt1, const cv::Point2f& pt2) {
    std::vector<double> lineParams;
    double a, b, c;
    _line_equation_by_points(pt1, pt2, a, b, c);
    lineParams.push_back(a);
    lineParams.push_back(b);
    lineParams.push_back(c);
    return lineParams;
}

static bool _are_lines_identical(double a1, double b1, double c1, double a2, double b2, double c2) {
    double a1b2 = a1*b2;
    double a2b1 = a2*b1;
    double a1c2 = a1*c2;
    double a2c1 = a2*c1;
    double b1c2 = b1*c2;
    double b2c1 = b2*c1;
    return (_almost_equal(a1b2, a2b1) && _almost_equal(a1c2, a2c1) && _almost_equal(b1c2, b2c1));
}

static bool _are_lines_identical(const std::vector<double>& side1Params, const std::vector<double>& side2Params, double sideCExtraParam) {
    return (_are_lines_identical(side1Params[0], side1Params[1], side1Params[2], side2Params[0], side2Params[1], side2Params[2]-sideCExtraParam) || _are_lines_identical(side1Params[0], side1Params[1], side1Params[2], side2Params[0], side2Params[1], side2Params[2]+sideCExtraParam));
}

static bool _line_intersect(double a1, double b1, double c1, double a2, double b2, double c2, cv::Point2f& intersectPt) {
    double det = a1*b2 - a2*b1;
    if (!_almost_equal(det, 0)) {
        intersectPt.x = static_cast<float>(-(b2*c1-b1*c2)/det);
        intersectPt.y = static_cast<float>(-(a1*c2-a2*c1)/det);
        return true;
    }
    return false;
}

static bool _line_intersect(const cv::Point2f& pt11, const cv::Point2f& pt12, const cv::Point2f& pt21, const cv::Point2f& pt22, cv::Point2f& intersectPt) {
    double a1 = pt12.y - pt11.y;
    double b1 = pt11.x - pt12.x;
    double c1 = pt12.x*pt11.y - pt11.x*pt12.y;
    
    double a2 = pt22.y - pt21.y;
    double b2 = pt21.x - pt22.x;
    double c2 = pt22.x*pt21.y - pt21.x*pt22.y;
    return _line_intersect(a1, b1, c1, a2, b2, c2, intersectPt);
}

static bool _are_intersect_lines(const std::vector<double>& side1Params, const std::vector<double>& side2Params, double sideCExtraParam, cv::Point2f& intersectPt1, cv::Point2f& intersectPt2) {
    return (_line_intersect(side1Params[0], side1Params[1], side1Params[2], side2Params[0], side2Params[1], side2Params[2]-sideCExtraParam, intersectPt1) && _line_intersect(side1Params[0], side1Params[1], side1Params[2], side2Params[0], side2Params[1], side2Params[2]+sideCExtraParam, intersectPt2));
}

static bool _find_gamma_intersect_points(const std::vector<cv::Point2f>& polygon, unsigned int c, unsigned int ptIdx, const cv::Point2f& side1StartVertex, const cv::Point2f& side1EndVertex, const cv::Point2f& side2StartVertex, const cv::Point2f& side2EndVertex, cv::Point2f& intersectPt1, cv::Point2f& intersectPt2) {
    std::vector<double> side1Params = _line_equation_parameters(side1StartVertex, side1EndVertex);
    std::vector<double> side2Params = _line_equation_parameters(side2StartVertex, side2EndVertex);
    
    double ptHeight = _height(ptIdx, polygon, c);
    double dist = sqrt(side2Params[0]*side2Params[0] + side2Params[1]*side2Params[1]);
    double sideCExtraParam = 2*ptHeight*dist;
    
    if (_are_intersect_lines(side1Params, side2Params, sideCExtraParam, intersectPt1, intersectPt2)) {
        return true;
    } else if (_are_lines_identical(side1Params, side2Params, sideCExtraParam)) {
        intersectPt1 = side1StartVertex;
        intersectPt2 = side1EndVertex;
        return true;
    }
    return false;
}

static int _sign(double num) {
    if (num > 0) {
        return 1;
    }
    if (num < 0) {
        return -1;
    }
    return 0;
}

static bool _are_on_the_same_side_of_line(const cv::Point2f& pt1, const cv::Point2f& pt2, const cv::Point2f& a, const cv::Point2f& b) {
    double a1, b1, c1;
    _line_equation_by_points(a, b, a1, b1, c1);
    double pt1OnLine = a1*pt1.x + b1*pt1.y + c1;
    double pt2OnLine = a1*pt2.x + b1*pt2.y + c1;
    return _sign(pt1OnLine) == _sign(pt2OnLine);
}

static bool _gamma(unsigned int ptIdx, cv::Point2f& gammaPt, const std::vector<cv::Point2f>& polygon, unsigned int a, unsigned int c) {
    unsigned int numOfPoints = polygon.size();
    cv::Point2f intersectPt1, intersectPt2;
    if (!_find_gamma_intersect_points(polygon, c, ptIdx, polygon[a], polygon[_predecessor(a, numOfPoints)], polygon[c], polygon[_predecessor(c, numOfPoints)], intersectPt1, intersectPt2)) {
        return false;
    }
    if (_are_on_the_same_side_of_line(intersectPt1, polygon[_successor(c, numOfPoints)], polygon[c], polygon[_predecessor(c, numOfPoints)])) {
        gammaPt = intersectPt1;
    } else {
        gammaPt = intersectPt2;
    }
    return true;
}

//////////////////////////////////////////////////////////////////////////////////////////////////////////////

#define INTERSECT_BELOW 1
#define INTERSECT_ABOVE 2
#define INTERSECT_CRITICAL 3

static double _angle_of_line_from_0x_axis(const cv::Point2f& a, const cv::Point2f& b) {
    double dy = b.y - a.y;
    double dx = b.x - a.x;
    double angle = std::atan2(dy, dx)*180/CV_PI;
    return (angle < 0) ? angle + 360 : angle;
}

static bool _less_or_equal(double num1, double num2) {
    return (num1 < num2 || _almost_equal(num1, num2));
}

static bool _is_angle_between(double angle1, double angle2, double angle3) {
    int angleDelta = (int)(angle2 - angle3);
    if (angleDelta %180 > 0) {
        return (angle3 < angle1 && angle1 < angle2);
    } else {
        return (angle2 < angle1 && angle1 < angle3);
    }
}

static bool _is_angle_between_non_reflex(double angle1, double angle2, double angle3) {
    if (std::abs(angle2 - angle3) > 180) {
        if (angle2 > angle3) {
            if (angle2 < angle1 && _less_or_equal(angle1, 360)) {
                return true;
            }
            if (_less_or_equal(0, angle1) && angle1 < angle3) {
                return true;
            }
            return false;
        } else {
            if (angle3 < angle1 && _less_or_equal(angle1, 360)) {
                return true;
            }
            if (_less_or_equal(0, angle1) && angle1 < angle2) {
                return true;
            }
            return false;
        }
    } else {
        return _is_angle_between(angle1, angle2, angle3);
    }
}

static double _opposite_angle(double angle) {
    if (angle > 180) {
        return angle - 180;
    }
    return angle + 180;
}

static bool _is_opposite_angle_between_non_reflex(double angle1, double angle2, double angle3) {
    double angle1Opposite = _opposite_angle(angle1);
    return _is_angle_between_non_reflex(angle1Opposite, angle2, angle3);
}

static bool _is_flush_angle_between_pred_and_succ(double& angleFlushEdge, double anglePred, double angleSucc) {
    if (_is_angle_between_non_reflex(angleFlushEdge, anglePred, angleSucc)) {
        return true;
    } else if (_is_opposite_angle_between_non_reflex(angleFlushEdge, anglePred, angleSucc)) {
        angleFlushEdge = _opposite_angle(angleFlushEdge);
        return true;
    }
    return false;
}

static bool _is_gamma_angle_between(double angleGamma, double angle1, double angle2) {
    return _is_angle_between_non_reflex(angleGamma, angle1, angle2);
}

static unsigned int _intersect_above_or_below(unsigned int succPredIdx, unsigned int ptIdx, const std::vector<cv::Point2f>& polygon, unsigned int c) {
    if (_height(succPredIdx, polygon, c) > _height(ptIdx, polygon, c)) {
        return INTERSECT_ABOVE;
    } else {
        return INTERSECT_BELOW;
    }
}

static bool _is_gamma_angle_equal_to(double gammaAngle, double angle) {
    return _almost_equal(gammaAngle, angle);
}

static unsigned int _intersect(double angleGammaAndPt, unsigned int ptIdx, const std::vector<cv::Point2f>& polygon, unsigned int c) {
    unsigned int numOfPoints = polygon.size();
    double anglePtPred = _angle_of_line_from_0x_axis(polygon[_predecessor(ptIdx, numOfPoints)], polygon[ptIdx]);
    double anglePtSucc = _angle_of_line_from_0x_axis(polygon[_successor(ptIdx, numOfPoints)], polygon[ptIdx]);
    double angleFlushEdge = _angle_of_line_from_0x_axis(polygon[_predecessor(c, numOfPoints)], polygon[c]);
    
    if (_is_flush_angle_between_pred_and_succ(angleFlushEdge, anglePtPred, anglePtSucc)) {
        if (_is_gamma_angle_between(angleGammaAndPt, anglePtPred, angleFlushEdge) || _almost_equal(angleGammaAndPt, anglePtPred)) {
            return _intersect_above_or_below(_predecessor(ptIdx, numOfPoints), ptIdx, polygon, c);
        } else if (_is_gamma_angle_between(angleGammaAndPt, anglePtSucc, angleFlushEdge) || _almost_equal(angleGammaAndPt, anglePtSucc)) {
            return _intersect_above_or_below(_successor(ptIdx, numOfPoints), ptIdx, polygon, c);
        }
    } else {
        if (_is_gamma_angle_between(angleGammaAndPt, anglePtPred, anglePtSucc) ||
            (_is_gamma_angle_equal_to(angleGammaAndPt, anglePtPred) && !_is_gamma_angle_equal_to(angleGammaAndPt, angleFlushEdge)) ||
            (_is_gamma_angle_equal_to(angleGammaAndPt, anglePtSucc) && !_is_gamma_angle_equal_to(angleGammaAndPt, angleFlushEdge))) {
            return INTERSECT_BELOW;
        }
    }
    return INTERSECT_CRITICAL;
}

static bool _intersect_below(const cv::Point2f& gammaPt, unsigned int ptIdx, const std::vector<cv::Point2f>& polygon, unsigned int c) {
    double angleOfGammaAndPt = _angle_of_line_from_0x_axis(polygon[ptIdx], gammaPt);
    return _intersect(angleOfGammaAndPt, ptIdx, polygon, c) == INTERSECT_BELOW;
}

static void _move_a_if_low_and_b_if_high(const std::vector<cv::Point2f>& polygon, unsigned int& a, unsigned int& b, unsigned int c) {
    unsigned int numOfPoints = polygon.size();
    cv::Point2f gammaA;
    while (_height(b, polygon, c) > _height(a, polygon, c)) {
        if (_gamma(a, gammaA, polygon, a, c) && _intersect_below(gammaA, b, polygon, c)) {
            _advance(b, numOfPoints);
        } else {
            _advance(a, numOfPoints);
        }
    }
}

//////////////////////////////////////////////////////////////////////////////////////////////////////////////

static void _search_for_b_tangency(const std::vector<cv::Point2f>& polygon, unsigned int a, unsigned int& b, unsigned int c) {
    unsigned int numOfPoints = polygon.size();
    cv::Point2f gammaOfB;
    while (_gamma(b, gammaOfB, polygon, a, c) && _intersect_below(gammaOfB, b, polygon, c) &&
        _greater_or_equal(_height(b, polygon, c), _height(_predecessor(a, numOfPoints), polygon, c))) {
        _advance(b, numOfPoints);
    }
}

static void _update_AC_sides(const std::vector<cv::Point2f>& polygon, unsigned int a, unsigned int c, cv::Point2f& sideAStartVertex, cv::Point2f& sideAEndVertex, cv::Point2f& sideCStartVertex, cv::Point2f& sideCEndVertex) {
    unsigned int numOfPoints = polygon.size();
    sideCStartVertex = polygon[_predecessor(c, numOfPoints)];
    sideCEndVertex = polygon[c];
    sideAStartVertex = polygon[_predecessor(a, numOfPoints)];
    sideAEndVertex = polygon[a];
}


//////////////////////////////////////////////////////////////////////////////////////////////////////////////

#define SIDE_A_TANGENT 0
#define SIDE_B_TANGENT 1
#define SIDES_FLUSH 2

static bool _intersects_above(const cv::Point2f& gammaPt, unsigned int ptIdx, const std::vector<cv::Point2f>& polygon, unsigned int c) {
    double angleOfGammaAndPt = _angle_of_line_from_0x_axis(gammaPt, polygon[ptIdx]);
    return _intersect(angleOfGammaAndPt, ptIdx, polygon, c) == INTERSECT_ABOVE;
}

static bool _is_not_b_tangency(const std::vector<cv::Point2f>& polygon, unsigned int a, unsigned int b, unsigned int c) {
    unsigned int numOfPoints = polygon.size();
    cv::Point2f gammaOfB;
    if ((_gamma(b, gammaOfB, polygon, a, c) && _intersects_above(gammaOfB, b, polygon, c)) || 
        (_height(b, polygon, c) < _height(_predecessor(a, numOfPoints), polygon, c))) {
        return true;
    }
    return false;
}

static cv::Point2f _middle_point(const cv::Point2f& pt1, const cv::Point2f& pt2) {
    double mx = static_cast<double>((pt1.x+pt2.x)/2.0);
    double my = static_cast<double>((pt1.y+pt2.y)/2.0);
    return cv::Point2f(static_cast<float>(mx), static_cast<float>(my));
}

static bool _middle_point_of_side_b(cv::Point2f& middlePointOfSideB, const cv::Point2f& sideAStartVertex, const cv::Point2f& sideAEndVertex, const cv::Point2f& sideBStartVertex, const cv::Point2f& sideBEndVertex, const cv::Point2f& sideCStartVertex, const cv::Point2f& sideCEndVertex) {
    cv::Point2f vertexA, vertexC;
    if (!_line_intersect(sideBStartVertex, sideBEndVertex, sideCStartVertex, sideCEndVertex, vertexA) ||
        !_line_intersect(sideBStartVertex, sideBEndVertex, sideAStartVertex, sideAEndVertex, vertexC)) {
        return false;
    }
    middlePointOfSideB = _middle_point(vertexA, vertexC);
    return true;
}

static cv::Point2f _find_vertex_c_on_side_b(const std::vector<cv::Point2f>& polygon, unsigned int a, unsigned int c, const cv::Point2f& sideBStartVertex,  cv::Point2f& sideBEndVertex, const cv::Point2f& sideCStartVertex, const cv::Point2f& sideCEndVertex) {
    unsigned int numOfPoints = polygon.size();
    cv::Point2f intersectPt1, intersectPt2;
    if (!_find_gamma_intersect_points(polygon, c, _predecessor(a, numOfPoints), sideBStartVertex, sideBEndVertex, sideCStartVertex, sideCEndVertex, intersectPt1, intersectPt2)) {
        // Error
        CV_Error(cv::Error::StsInternal, "Error Vertex C on Side B");
    }
    if (_are_on_the_same_side_of_line(intersectPt1, polygon[_successor(c, numOfPoints)], polygon[c], polygon[_predecessor(c, numOfPoints)])) {
        return intersectPt1;
    } else {
        return intersectPt2;
    }
}

static void _update_AB_sides(const std::vector<cv::Point2f>& polygon, unsigned int a, unsigned int b, int c, unsigned int& validationFlag, cv::Point2f& sideAStartVertex, cv::Point2f& sideAEndVertex, cv::Point2f& sideBStartVertex, cv::Point2f& sideBEndVertex, const cv::Point2f& sideCStartVertex, const cv::Point2f& sideCEndVertex) {
    unsigned int numOfPoints = polygon.size();
    sideBStartVertex = polygon[_predecessor(b, numOfPoints)];
    sideBEndVertex = polygon[b];
    
    cv::Point2f sideBMiddlePoint;
    if (_middle_point_of_side_b(sideBMiddlePoint, sideAStartVertex, sideAEndVertex, sideBStartVertex, sideBEndVertex, sideCStartVertex, sideCEndVertex) && (_height(sideBMiddlePoint, polygon, c) < _height(_predecessor(a, numOfPoints), polygon, c))) {
        sideAStartVertex = polygon[_predecessor(a, numOfPoints)];
        sideAEndVertex = _find_vertex_c_on_side_b(polygon, a, c, sideBStartVertex, sideBEndVertex, sideCStartVertex, sideCEndVertex);
        validationFlag = SIDE_A_TANGENT;
    } else {
        validationFlag = SIDES_FLUSH;
    }
}

static void _update_b_side(const std::vector<cv::Point2f>& polygon, unsigned int a, unsigned int b, unsigned int c, unsigned int& validationFlag, cv::Point2f& sideBStartVertex, cv::Point2f& sideBEndVertex) {
    if (!_gamma(b, sideBStartVertex, polygon, a, c)) {
        CV_Error(cv::Error::StsInternal, "Error Side B Gamma");
    }
    sideBEndVertex = polygon[b];
    validationFlag = SIDE_B_TANGENT;
}
                                            
//////////////////////////////////////////////////////////////////////////////////////////////////////////////

static double _distance_between_points(const cv::Point2f& pt1, const cv::Point2f& pt2) {
    double dx = pt1.x - pt2.x;
    double dy = pt1.y - pt2.y;
    return std::sqrt(dx*dx+dy*dy);
}

static bool _is_point_on_line_segment(const cv::Point2f& pt, const cv::Point2f& linePt1, const cv::Point2f& linePt2) {
    double d1 = _distance_between_points(pt, linePt1);
    double d2 = _distance_between_points(pt, linePt2);
    double lineLen = _distance_between_points(linePt1, linePt2);
    return _almost_equal(d1+d2, lineLen);
}

static bool _is_valid_min_triangle(const cv::Point2f& vertexA, const cv::Point2f& vertexB, const cv::Point2f& vertexC, const std::vector<cv::Point2f>& polygon, unsigned int a, unsigned int b, unsigned int validationFlag, const cv::Point2f& sideAStartVertex, const cv::Point2f& sideAEndVertex, const cv::Point2f& sideBStartVertex, const cv::Point2f& sideBEndVertex, const cv::Point2f& sideCStartVertex, const cv::Point2f& sideCEndVertex) {
    cv::Point2f middlePtSideA = _middle_point(vertexB, vertexC);
    cv::Point2f middlePtSideB = _middle_point(vertexA, vertexC);
    cv::Point2f middlePtSideC = _middle_point(vertexA, vertexB);
    
    unsigned int numOfPoints = polygon.size();
    bool sideAValid = (validationFlag == SIDE_A_TANGENT ? 
        _are_equal_points(middlePtSideA, polygon[_predecessor(a, numOfPoints)]) :
        _is_point_on_line_segment(middlePtSideA, sideAStartVertex, sideAEndVertex)
    );
    bool sideBValid = (validationFlag == SIDE_B_TANGENT ? 
        _are_equal_points(middlePtSideB, polygon[b]) :
        _is_point_on_line_segment(middlePtSideB, sideBStartVertex, sideBEndVertex)
    );
    bool sideCValid = (validationFlag == SIDES_FLUSH || _is_point_on_line_segment(middlePtSideC, sideCStartVertex, sideCEndVertex));
    
    return (sideAValid && sideBValid && sideCValid);
}

static bool _is_local_min_triangle(cv::Point2f& vertexA, cv::Point2f& vertexB, cv::Point2f& vertexC, const std::vector<cv::Point2f>& polygon, unsigned int a, unsigned int b, unsigned int validationFlag, const cv::Point2f& sideAStartVertex, const cv::Point2f& sideAEndVertex, const cv::Point2f& sideBStartVertex, const cv::Point2f& sideBEndVertex, const cv::Point2f& sideCStartVertex, const cv::Point2f& sideCEndVertex) {
    if (!_line_intersect(sideAStartVertex, sideAEndVertex, sideBStartVertex, sideBEndVertex, vertexC) ||
        !_line_intersect(sideAStartVertex, sideAEndVertex, sideCStartVertex, sideCEndVertex, vertexB) ||
        !_line_intersect(sideBStartVertex, sideBEndVertex, sideCStartVertex, sideCEndVertex, vertexA)) {
        return false;
    }
    return _is_valid_min_triangle(vertexA, vertexB, vertexC, polygon, a, b, validationFlag, sideAStartVertex, sideAEndVertex, sideBStartVertex, sideBEndVertex, sideCStartVertex, sideCEndVertex);
}

static double _area_of_triangle(const cv::Point2f& p1, const cv::Point2f& p2, const cv::Point2f& p3) {
    double posTerm = p1.x*p2.y + p1.y*p3.x + p2.x*p3.y;
    double negTerm = p2.y*p3.x + p1.x*p3.y + p1.y*p2.x;
    double det = posTerm - negTerm;
    return std::abs(det)/2;
}

static void _update_min_area_enclosing_triangle(std::vector<cv::Point2f>& triangle, double& area, const cv::Point2f& vertexA, const cv::Point2f& vertexB, const cv::Point2f& vertexC) {
    double triArea = _area_of_triangle(vertexA, vertexB, vertexC);
    if (triArea < area) {
        triangle.clear();
        triangle.push_back(vertexA);
        triangle.push_back(vertexB);
        triangle.push_back(vertexC);
        area = triArea;
    }
}

//////////////////////////////////////////////////////////////////////////////////////////////////////////////

static void _find_min_enclosing_triangle(const std::vector<cv::Point2f>& points, std::vector<cv::Point2f>& triangle, double& area) {
    std::vector<cv::Point2f> polygon;
    bool clockwise = true;
    bool returnPoints = true;
    cv::convexHull(points, polygon, clockwise, returnPoints);
    
    area = std::numeric_limits<double>::max();
    triangle.clear();
    
    unsigned int numOfPoints = polygon.size();
    if (numOfPoints == 0) {
        return;
    } else if (numOfPoints <= 3) {
        if (numOfPoints == 1) {
            triangle.push_back(polygon[0]);
            triangle.push_back(polygon[0]);
            triangle.push_back(polygon[0]);
        } else if (numOfPoints == 2) {
            triangle.push_back(polygon[0]);
            triangle.push_back(polygon[1]);
            triangle.push_back(polygon[0]);
        } else {
            triangle.push_back(polygon[0]);
            triangle.push_back(polygon[1]);
            triangle.push_back(polygon[2]);
        }
        area = _area_of_triangle(triangle[0], triangle[1], triangle[2]);
    } else {
        cv::Point2f vertexA, vertexB, vertexC;
        cv::Point2f sideAStartVertex, sideAEndVertex, sideBStartVertex, sideBEndVertex, sideCStartVertex, sideCEndVertex;
        
        unsigned int a = 1, b = 2, c, validationFlag;
        for (c = 0; c < numOfPoints; c++) {
            _advance_b_to_right_chain(polygon, b, c);
            _move_a_if_low_and_b_if_high(polygon, a, b, c);
            _search_for_b_tangency(polygon, a, b, c);
            _update_AC_sides(polygon, a, c, sideAStartVertex, sideAEndVertex, sideCStartVertex, sideCEndVertex);
            
            if (_is_not_b_tangency(polygon, a, b, c)) {
                _update_AB_sides(polygon, a, b, c, validationFlag, sideAStartVertex, sideAEndVertex, sideBStartVertex, sideBEndVertex, sideCStartVertex, sideCEndVertex);
            } else {
                _update_b_side(polygon, a, b, c, validationFlag, sideBStartVertex, sideBEndVertex);
            }
            
            if (_is_local_min_triangle(vertexA, vertexB, vertexC, polygon, a, b, validationFlag, sideAStartVertex, sideAEndVertex, sideBStartVertex, sideBEndVertex, sideCStartVertex, sideCEndVertex)) {
                _update_min_area_enclosing_triangle(triangle, area, vertexA, vertexB, vertexC);
            }
        }
    }
}


int main(int argc, char **argv) {
    cv::Mat src = cv::imread("/home/xlll/Downloads/opencv/samples/data/lena.jpg", cv::IMREAD_COLOR);
    if (src.empty()) {
        std::cout << "failed to read image!" << std::endl;
        return EXIT_FAILURE;
    }
    
    cv::Mat binary;
    cv::Canny(src, binary, 250, 1000);
    
    std::vector<std::vector<cv::Point>> contours;
    cv::findContours(binary, contours, cv::RETR_LIST, cv::CHAIN_APPROX_NONE);
    std::vector<std::vector<cv::Point2f>> triangles, triangles2;
    cv::Mat src2 = src.clone();
    for (int i = 0; i < contours.size(); i++) {
        /*
         * std::vector<cv::Point2f> triangle;
        cv::minEnclosingTriangle(contours[i], triangle);
        triangles.push_back(triangle);
        for (int j = 0; j < contours[i].size(); j++) {
            cv::circle(src, contours[i][j], 1, cv::Scalar(255, 0, 0), cv::FILLED);
        }
        */
        
        std::vector<cv::Point2f> points;
        for (int j = 0; j < contours[i].size(); j++) {
            points.push_back(cv::Point2f(static_cast<float>(contours[i][j].x), static_cast<float>(contours[i][j].y)));
        }
        std::vector<cv::Point2f> triangle, triangle2;
        cv::minEnclosingTriangle(points, triangle);
        triangles.push_back(triangle);
        for (int j = 0; j < points.size(); j++) {
            cv::circle(src, points[j], 1, cv::Scalar(255, 0, 0), cv::FILLED);
            cv::circle(src2, points[j], 1, cv::Scalar(255, 0, 0), cv::FILLED);
        }
        double area;
        _find_min_enclosing_triangle(points, triangle2, area);
        triangles2.push_back(triangle2);
    }
    for (int i = 0; i < triangles.size(); i++) {
        for (int j = 0; j < triangles[i].size(); j++) {
            cv::line(src, triangles[i][j], triangles[i][(j+1)%3], cv::Scalar(0, 255, 0), 1, cv::LINE_AA);
            cv::line(src2, triangles2[i][j], triangles2[i][(j+1)%3], cv::Scalar(0, 255, 0), 1, cv::LINE_AA);
            
            std::cout << triangles[i][j] << " -- " << triangles2[i][j] << std::endl;
        }
    }
    
    cv::imshow("src", src);
    cv::imshow("src2", src2);
    cv::imshow("binary", binary);
    cv::waitKey(0);
    
    return 0;
}
